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The reaction-diffusion master equation (RDME) is a lattice stochastic reaction-diffusion model 
that has been used to study spatially distributed cellular processes. The RDME has been shown to 
have the drawback of losing bimolecular reactions in the continuum limit that the lattice spacing 
approaches zero (in two or more dimensions). In this work we derive a new convergent RDME 
(CRDME) that eliminates this problem. The CRDME is obtained by finite volume discretization of 
a spatially-continuous stochastic reaction-diffusion model. We demonstrate the empirical numerical 
convergence of reaction time statistics associated with the CRDME. Although the reaction time 
statistics of the RDME diverge as the lattice spacing approaches zero, we show they approach those 
of the CRDME for sufficiently large lattice spacings or slow bimolecular reaction rates. As such, 
the RDME may be interpreted as an approximation to the CRDME in several asymptotic limits. 



I. INTRODUCTION 

Computational models of biochemical systems within 
individual cells have become a common tool used in 
studying cellular processes and behavior J55J 1311 I3H E2] • 
The spatially distributed nature of chemical pathways 
inside cells may, in certain cases, be more accurately 
modeled by including the explicit spatial movement of 
molecules. Examples of such processes include whether 
signals can propagate from the plasma membrane to nu- 
cleus |33j , how cell shape can modify information flow in 
signaling networks [34 , how the variable density of chro- 
matin influences the search time of proteins for DNA 
binding sites [55], or how different regions of cytosolic 
space may separate to different chemical states [117] . 

One important consideration in developing these mod- 
els is that at the scale of a single cell many biochemical 
processes are stochastic [5J 0] [35] . Stochastic reaction- 
diffusion models have been used to account for the 
stochasticity inherent in the chemical reaction process 
and the diffusion of proteins and mRNAs. These mod- 
els approximate individual molecules as points diffusing 
within cells. They are more macroscopic descriptions 
than quantum mechanical or molecular dynamics mod- 
els, which can resolve detailed interactions between a few 
molecules on timescales of milliseconds [35]. They are 
more microscopic descriptions than deterministic three- 
dimensional reaction-diffusion PDEs for the average con- 
centration of each species of molecule. 

Three stochastic reaction-diffusion models that have 
been used to study cellular processes are: the Doi 
model [51 [7J |32] , the Smoluchowski diffusion limited re- 
action model (SDLR) [211 [37] , and the reaction-diffusion 
master equation (RDME) [Efl EH H3 EH E3 ED| • In the 
Doi model [6j [Jj positions of molecules are represented 
as points undergoing Brownian motion. Bimolecular re- 
actions between two molecules occur with a fixed proba- 
bility per unit time when two reactants are separated by 
less than some specified "reaction radius". The SDLR 
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model differs by representing bimolecular reactions in one 
of two ways; either occurring instantaneously, or with 
fixed probability per unit time, when two reactants' sep- 
aration is exactly the reaction-radius [291 I37j . In each 
of the models unimolccular reactions represent internal 
processes. They are assumed to occur with exponen- 
tially distributed times based on a specified reaction-rate 
constant. 

The RDME can be interpreted as an extension of the 
non-spatial chemical master equation (CME) [151 [T71 1351 
HO] model for stochastic chemical kinetics. In the RDME 
space is partitioned by a mesh into a collection of vox- 
els. The diffusion of molecules is modeled as a continuous 
time random walk on the mesh, with bimolecular reac- 
tions occurring with a fixed probability per unit time 
for molecules within the same voxel. Within each voxel 
molecules are assumed well- mixed (i.e. uniformly dis- 
tributed), and molecules of the same species are indis- 
tinguishable. Mathematically, the RDME is the forward 
Kolmogorov equation for a continuous-time jump Markov 
process. 

In each of the three models the state of the chemi- 
cal system is given by stochastic processes for the num- 
ber of molecules of each chemical species and the corre- 
sponding positions of each molecule. There are two dif- 
ferent mathematical formulations of each model; either 
coupled systems of equations for the probability densi- 
ties of having a given state at a specified time, or equa- 
tions for the evolution of the stochastic processes them- 
selves. The former description leads to, possibly infinite, 
coupled systems of partial integral differential equations 
(PIDEs) for the Doi/SDLR models, and, possibly infi- 
nite, coupled systems of ordinary differential equations 
(ODEs) for the RDME. The high dimensionality of these 
equations for typical biological systems prevents the use 
of standard numerical methods for PDEs/ODEs in low- 
dimensions. Instead, the probability density solutions 
to these equations are approximated by simulating the 
underlying stochastic processes by Monte Carlo meth- 
ods pa Bang nana Eg. 

While the RDME can be used as an independent 
model, in applications it is often interpreted as a for- 
mal approximation of the Doi or SDLR models [13]. In 
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practice, using the RDME to approximate either of the 
Doi or SDLR models can be problematic as it has been 
shown that in the continuum limit where the mesh spac- 
ing in the RDME is taken to zero bimolecular reactions 
are lost [3TJI2I] (in two or more dimensions). That is, the 
time at which two molecules will react becomes infinite as 
the mesh spacing is taken to zero. In general the RDME 
can only provide a reasonable approximation to the Doi 
or SDLR models for mesh spacings that are neither too 
large or too small [25]. This error in this approximation 
can not be made arbitrarily small [25] . 

In this work we discretize the Doi model so as to obtain 
a forward Kolmogorov equation describing a continuous- 
time jump Markov process similar to the RDME. We seek 
a discretization of this form so that we may use Monte 
Carlo methods to simulate the underlying stochastic pro- 
cess for systems of many reacting molecules. Our con- 
vergent RDME (CRDME) approximates the Doi model, 
but retains bimolecular reactions as the mesh spacing ap- 
proaches zero. In particular, for T the random variable 
for the time at which two molecules react, we show in 
two-dimensions by numerical simulation that both the 
survival time distribution, Pr [T > t] , and the mean re- 
action time converge to finite values as the mesh spac- 
ing approaches zero in the CRDME (in contrast to the 
RDME). 

The new CRDME retains many of the benefits of 
the original RDME model, and allows the re-use of the 
many extensions of the RDME that have been devel- 
oped. Examples of these extensions include Cartesian- 
grid methods for complex geometries |27j ; the incorpora- 
tion of drift due to potentials [251 H3j ; advective velocity 
fields [20 ; non-Cartesian meshes [TT]; GPU optimized 
simulation methods |41j ; adaptive mesh refinement tech- 
niques (AMR) to improve the approximation of molecu- 
lar diffusion [3J ; and multiscale couplings to more macro- 
scopic models [15) . In addition to retaining bimolecular 
reactions as the mesh spacing approaches zero, we show 
that the CRDME is approximable by the RDME for ap- 
propriate mesh spacings and parameter choices. 



species A when the total number of molecules of species 
A is a. The state vector of the species A molecules is 
then given by q a = (qf , . . . , q a a ) € R da . Define q b and q c 
similarly. We denote by f^ a ' b ^{q a ,q b ,q c ^t) the proba- 
bility density for there to be a molecules of species A, b 
molecules of species B, and c molecules of species C at 
time t located at the positions q a , q b , and q c . Molecules 
of the same species are assumed indistinguishable. The 
evolution of /( a ' fc > c ' is given by 
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dt 



(q a , q b : q c , t) = (L + R) (q a ,q b , q c , t) 



(1) 

Note, with the subsequent definitions of the operators L 
and R this will give a coupled system of PIDEs over all 
possible values of (a, b, c). More general systems that al- 
low unbounded production of certain species would result 
in an infinite number of coupled PIDEs. The diffusion 
operator, L, is defined by 
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where A" denotes the Laplacian in the coordinate qf and 
D A the diffusion constant of species A. D B , D c , A^, and 
A^j are defined similarly. 

To define the reaction operator, R, we introduce no- 
tations for removing or adding a specific molecule to the 
state q a . Let 



q a \q? = {q a i, 

q a yjq= («?,. 



■ ,qU,qi +1 ,...,ql) 



q a \ q will denote q a with any one component with 
the value q removed. Denote by lro lT . b i(r) the indica- 
tor function of the interval [0, n,], and let Bf = {q e 
\q ~ qf\ < r b/2} label the set of points a reactant 
could be at to produce a molecule of species C at qf. 
The Doi reaction operator, R, is then 
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II. GENERAL DOI AND RDME MODELS 

We first illustrate how the bimolecular reaction A + 
B — > C would be described by the Doi model and the 
standard RDME in K d . In the Doi model, bimolecular 
reactions are characterized by two parameters; the sepa- 
ration at which molecules may begin to react, n,, and the 
probability per unit time the molecules react when within 
this separation, A. When a molecule of species A and a 
molecule of species B react we assume the C molecule 
they produce is placed midway between them. Note the 
important point that in each of the models molecules are 
modeled as points. 

We now formulate the Doi model as an infinite coupled 
system of partial integral differential equations (PIDEs) . 
Let qf g W 1 denote the position of the Ith molecule of 
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We now describe the RDME, in a form we derived 
in |23] that has the advantage of representing a chemi- 
cal system's state in a similar manner to the Doi model. 
Using this representation allows for easier comparison of 
the RDME and Doi models. Partition R d into a Carte- 
sian lattice of voxels with width h and hypervolume h d . 
When in the same voxel, an A and B molecule may react 
with probability per unit time k/h d . Here k represents 
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the macroscopic bimolecular reaction-rate constant of the 
reaction A + B — > C, with units of hypervolume per unit 
time. Let jf £ 1 d denote the multi-index of the voxel 
centered at hjf that contains the Zth molecule of species 
A when there are a molecules of species A. The position 
of the molecule is assumed to be uniformly distributed 
(i.e. well-mixed) within this voxel. Let j a — . . . , j£) 
denote the state vector for the voxels containing the a 
molecules of species A. Define j b and j c similarly, and 
let F{ i a ' b ' c \j a ,j b 1 j c ,t) denote the probability that there 
are (a, b, c) molecules of species A, B, and C at time t in 
the voxels given by j a , j b , and f. The RDME is then 
the coupled system of ODEs over all possible values for 
a, b, c, j a , j b , and j c , 



' /F;; """' ) (/'././•/) = (L h + R h )Ft' b ' c) (.f ././•/) , 

(4) 

where is a discretized approximation to L given by 

L h Ft b ' c) (3 a ,j b ,f,t) = 

(D A Al + D^A b h + D c Al)F^ M (f,j b J c ,t) . 

Here denotes the standard da-dimensional discrete 
Laplacian acting in the j a coordinate. We define the 
"standard" d-dimensional discrete Laplacian acting on a 
mesh function, f(j) on Z d , by 

d 

fc=i 

where denotes a unit vector along the kth coordinate 
axis of R d . 

The reaction operator, Rh, is given by 
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where ^h(if — il) denotes the Kronecker delta function 
equal to one when jf = j m and zero otherwise. Note 
that these equations are a formal discrete approximation 
to the Doi model, where two molecules may now react 
when in the same voxel with rate k/h d . 

We have shown that the RDME Q may be interpreted 
as a formal approximation to both Doi-like and Smolu- 
chowski models [53H2S] . As we proved in [53] and studied 
numerically in [2S], the RDME loses bimolecular reac- 
tions as h — > 0. This was shown rigorously for d = 3, 
where the divergence was shown to be like h . A sim- 
ple modification of our argument shows bimolecular reac- 
tions are lost for all d > 1, with a divergence like ln(h) for 



d = 2, and like hr d+2 for d > 2. More recently asymp- 
totic expansions were used in [21] to show the mean re- 
action time becomes infinite with the preceding rates as 
h — > (for d = 2 or d = 3). This loss of reaction occurs 
because molecules are modeled by points, and as h — > 
each voxel of the mesh shrinks to a point. Since two 
molecules must be in the same voxel to react, and in two 
or more dimensions two points can not find each other 
by diffusion, bimolecular reactions will never occur. 

It should be noted that while the RDME loses bimolec- 
ular reactions in the limit that h — > 0, we have shown that 
for fixed values of h it gives an asymptotic approximation 
in the binding radius, fb, to the SDLR model (2U |2"S"] . 
How accurate this approximation can be made is de- 
pendent on domain geometry and the parameters of the 
underlying chemical system [25 . In particular, h must 
be chosen sufficiently large that reactions within a voxel 
can be approximated by a well-mixed reaction with rate 
k/h d , while chosen sufficiently small that the diffusion of 
the molecules is well-approximated by a continuous time 
random walk on the mesh [I0l [24] . 

In the next section we develop a new convergent 
RDME (CRDME) to overcome these limitations of the 
standard RDME @. 



III. A CONVERGENT RDME (CRDME) 

To construct a convergent RDME (CRDME) we use 
a finite volume discretization of the Doi PIDEs ([!]) . For 
brevity we illustrate our approach on a simplified ver- 
sion of when there is only one molecule of A and one 
molecule of B in the system which may undergo the anni- 
hilation reaction A|B-> 0. The approach we describe 
can be extended to general multi-particle systems as bi- 
molecular reactions in the Doi model only involve multi- 
ple two-particle interactions of the same form, see (fTl). 

For now we work in d-dimensional free-space, K (for 
most biological models d = 2 or d = 3). Denote by 
x G R d the position of the molecule of species A and 
by y £ R d the position of the molecule of species B. 
In the Doi model these molecules diffuse independently, 
and may react with probability per unit time A when 
within a separation r^. (r^ is usually called the reaction- 
radius.) We let 1Z = {(x,y) \\x — y\< r^}, and denote 
the indicator function of this set by l^j. ( | a; — y\). The 
diffusion constants of the two molecules will be given by 
D A and D B respectively. 

Finally, we denote by p(x, y, t) the probability density 
the two molecules have not reacted and are at the posi- 
tions x and y at time t. Then reduces to 



dp 
di 



(x, y, t) = {D A A X + D B A y )p{x 7 y, t) 



-\tn(\x-y\)p(x,y,t). (7) 



(Here we have dropped the equation for the state a = 0, 
b = 0.) 
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We now show how to construct a new type of RDME 
by discretization of this equation. Note, while ^ can be 
solved analytically by switching to the separation coor- 
dinate, x — y, such approaches will not work for more 
general chemical systems, such as ([T]). As such, we illus- 
trate our CRDME ideas on ([7]), though ultimately they 
are intended for use in modeling more complex chemical 
systems with many chemical species, reactions, and com- 
plex geometries. For % E Z d and j E 7L d we partition M. 2d 
into a Cartesian grid of hypercubes, labeled by Vij. We 
assume Vij has coordinate-axis aligned edges with length 
h and center (cci,y_-) = (ih 7 jh). We denote the hyper- 
volume of a set, S 1 , by |5|. For example, the hypervolumc 
of the hypercube Vij is \Vij \ — h 2d . 

We make the approximation that the probability 
(x,y) E Vij at time t is given by Pij(t) = p(xi,yj,t) x 
| Vij | . Note that Pij (t) equivalently gives the probability 
that the particle of species A is in the ci-dimensional hy- 
percube of length h centered about ih, while the particle 
of species B is in the d-dimensional hypercube of length h 
about the point jh. Using this assumption we construct 
a finite volume discretization of ([7]) by integrating both 
sides of ([7]) over the hypercube Vij. As we did in [27] , 
we make the standard finite volume approximations for 
the integrals involving A x p and A y p to obtain discrete 
Laplacians given by ^ in the x and y coordinates. The 
reaction term is approximated by 

X l n {\x-y\)p{x,y,t)dxdy 

JViA 

^—Pi.jii) [ t n (\x-y\)dxdy 

y ij I J V i:i 
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Define the volume fraction 
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discretization then represents a new RDME for the two- 
particle system, given by the coupled system of ODEs 
over all values of E l? d , 



dP 



dt 



(t) = L h Pij(t) - XfaPijit). 



(8) 



Here L h = {D A L A + D B L B ), with L A and L? denoting 
standard d-dimensional discrete Laplacians (pi) (in the i 
and j coordinates respectively). By choosing an appro- 
priate discretization we have obtained an equation that 
has the form of the forward Kolomogorov equation for a 
continuous-time jump Markov process. We may therefore 
interpret the coefficients in ^ as transition rates, also 
called propensities, between states of the stochas- 

tic process. We subsequently refer to this equation as 
the CRDME. In the CRDME diffusion is handled in ex- 
actly the same manner as for the RDME. In contrast, the 
reaction mechanism in ^ now allows molecules in dis- 
tinct voxels, i and j, to react with a potentially non-zero 
probability per unit time, Xtfiij. 

As discussed earlier, the RDME is only physically valid 
when h is chosen sufficiently large that the timescale for 



two molecules to become uniformly distributed within a 
voxel by diffusion is much faster than that for a well- 
mixed bimolecular reaction to occur between them. By 
allowing molecules to react when in nearby voxels, our 
CRDME provides a correction when h is sufficiently small 
that this condition is violated. For > h molecules may 
potentially react when separated by multiple voxels, with 
the number of voxels apart two molecules can be and still 
react increasing as h — > 0. 



IV. EMPIRICAL CONVERGENCE OF THE 
CRDME 

We now demonstrate that in two-dimensions (d = 2) 
the survival time distribution and the mean reaction time 
for the CRDME ([8| converge to finite values as h — > 0. 
We show that in the corresponding RDME model the sur- 
vival time and mean reaction time diverge to 00 as h — > 0. 
In contrast, in the opposite limit that r^/h — > we 
demonstrate that the mean reaction time in the RDME 
approaches that of the CRDME. As is a PDE in four- 
dimensions, we do not directly solve the corresponding 
system of ODEs given by the CRDME @. Instead, we 
simulate the corresponding stochastic jump process for 
the molecule's motion and reaction by the well-known 
exact stochastic simulation algorithm (SSA) (also known 
as the Gillespie method [18 or kinetic Monte Carlo [S]). 

We assume each molecule moves within a square with 
sides of length L, ft — [0, L] x [0, L], with zero Neumann 
boundary conditions. These boundary conditions are en- 
forced by setting the transition rate for a molecule to 
hop from a given mesh voxel outside the domain to zero. 
For a specified number of mesh voxels, N, we discretize 
O into a Cartesian grid of squares with sides of length 
h = L/N. We assume D A — D B = D. Unless otherwise 
specified, all spatial units will be micrometers, with time 
in units of seconds. The SSA-based simulation algorithm 
can be summarized as: 

1. Specify D, L, N, A, and as input. 

2. Calculate <f>Qj. (See the supplementary material.) 

3. We assume the molecules are well-mixed at time 
t = 0. That is, the initial position of each molecule 
is sampled from a uniform distribution among all 
voxels of the mesh. 

4. Sample a time and direction of the next spatial hop 
by one of the molecules. (From ^ each molecule 
may hop to a neighbor in the x or y direction with 
probability per unit time D/h 2 .) 

5. Assuming the molecules are at i and j, if tf>ij =/= 
sample the next reaction time using the transition 
rate \4>%j- 

6. Select the smaller of the hopping and reaction 
times, and execute that event. Update the current 
time to the time of the event. 

7. If a reaction occurs the simulation ends. If a spatial 
hop occurs, return to 4. 

For all simulations we chose L = .2 /im. With this 
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FIG. 1: Survival time distributions vs. t for the CRDME and RDME when A = 10 9 s _1 . Each curve was estimated 
from 128000 simulations. The legends give the ratio, r^/h, used for each curve (note that = lnm was fixed and h 
successively halved). We see the convergence of the survival time distributions for the CRDME (up to sampling 
error), while the survival time distributions in the RDME diverge. 




FIG. 2: Mean reaction times vs. r^/h as h decreases by factors of two. Legends give the value of A for each curve 
(with units of s _1 ). Each mean reaction time was estimated from 128000 simulations. Note that 95% confidence 
intervals are drawn on each data point. (For some points they are smaller than the marker labeling the point.) Since 
we use a logarithmic x-axis, we see that for h sufficiently small the mean reaction time in the standard RDME 
diverges like ln(/i), while in the CRDME the mean reaction time converges to a finite value. 



choice the domain could be interpreted as small patch 
of membrane within a cell. A diffusion constant of 
D = 10/xm 2 s _1 was used for each molecule. The re- 
action radius, rt>, was chosen to be lnm. While physi- 
cal reaction radii are generally not measured experimen- 
tally, this choice falls between the measured width of the 
LexA DNA binding potential (« 5 angstroms [3D]) and 
the 5 nm reaction radius used for interacting membrane 
proteins in [jj]. 



We also simulated the stochastic process described by 
the corresponding RDME model. The bimolecular reac- 
tion rate was chosen to be k — Xirr^ to illustrate how 
the RDME approximates the CRDME as r h /h 0, but 
diverges as r^/h — > oo. Our motivation for this choice is 
explained in the next section. The simulation algorithm 
was identical to that just described, except that step 2 
was removed and step 5 modified so that two molecules 
could only react when within the same voxel (with prob- 
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ability per unit time k — \irr 2 /h 2 
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FIG. 3: Difference between successive points on the 
A = 10 9 CRDME curve in Fig. 2 vs. r h /h. The smaller 
of the two h values is used to label each point. The first 
and second order curves scale like h and h 2 respectively. 

Observe that the effective convergence rate to zero is 
closer to 0(h 2 ) than 0(h). 

Let T denote the random variable for the time at which 
the two molecules react. The survival time distribution, 
Pr [T >t], is 

Pr[T>t] = P {x,y,t)dxdy. 
Jo Jo 

Note that the reaction time distribution, Pr [T < t] — 
1 — Pr [T > t] . We estimate Pr [T > t] from the numer- 
ically sampled reaction times using the MATLAB ecdf 
function. In Fig. 1 we show the convergence (to within 
sampling error) of the estimated survival time distri- 
butions for the CRDME (right figure) as h — > (for 
A = 10 9 s _1 ). In the left figure we show the divergence 
as h — > of the estimated survival time distributions in 
the RDME. The continuing rightward shift of the dis- 
tribution as the mesh width is decreased to twenty times 
finer than the reaction radius shows the divergence of the 
reaction time to infinity. 

The mean reaction time, E[T], is given by 



E[T] 



Pr [T > t] dt. 



We estimated E[T] from the numerically sampled reac- 
tion times by calculating the sample mean. Fig. 2 shows 
the estimated mean reaction times for the CRDME and 
RDME models as A and rh/h are varied (note the cc-axis 
is logarithmic in rh/h). We see that as h — > the sam- 
pled mean reaction times in the CRDME converge to a 
fixed value. The rate of convergence in the CRDME for 
A = 10 9 is illustrated in Fig. 3. There we plot the differ- 
ence between successive estimated mean reaction times 
as h is halved. For small values of h this difference is 



seen to converge close to second order (as illustrated by 
the slope of the solid blue line) . 

In contrast, Fig. 2 shows that the sampled mean reac- 
tion time in the RDME diverges like ln(h) as discussed 
in [STJ [21] . For all A values the sampled mean reaction 
time in the RDME converges to that of the CRDME as 
rb/h — > 0. As A is decreased we see agreement between 
the RDME and CRDME for a larger range oir^/h values. 

Figs. 1, 2, and 3 demonstrate that, in contrast to the 
RDME, the reaction time statistics in the CRDME con- 
verge as h — > 0. That said, as rh/h — > we see that 
reaction time statistics of the RDME converge to those 
in the CRDME. Hence we may interpret the RDME as an 
approximation to the CRDME for rh/h <C 1. The accu- 
racy of using this approximation to describe the reaction- 
diffusion process will then depend on the relative sizes of 
A, D, and as we discuss in the next section (and illus- 
trated in Fig. 2). 



V. RDME AS AN APPROXIMATION OF THE 
CRDME FOR r b /h < 1 

We now show that the RDME may be interpreted 
as an asymptotic approximation to the CRDME for 
rh/h <C 1, with the accuracy of this approximation de- 
pending on the size of D, 7"b, and A. In the standard 
RDME two molecules can only react when within the 
same d-dimcnsional voxel. (If k denotes a macroscopic 
bimolecular reaction rate, the probability per unit time 
the molecules react is usually chosen to be k/h d , see §6§.) 
In contrast, our new model allows two molecules to re- 
act when in neighboring voxels. Even for large values of 
h, the volume fraction tfiij will be non-zero when i and 
j are neighboring voxels. That said, for j ^ i the vol- 
ume fraction </>jj will approach zero quicker as r^/h — > 
than (f>a. This relationship is shown in two-dimensions 
(d = 2) in Fig. 4. (We describe how the area fractions 
were calculated in the supplementary material.) 

We therefore expect that, asymptotically, when 
rh/h — > the particles will effectively only react when 
j = i. In this case 



\HnVu\ = I dydx 

I dy dx 

l[-h/2,h/2] d J{y\\ X -y\<r b } 
= \Br b \h d , 

where \B rh \ denotes the volume of the ei-dimensional 
sphere of radius tv The reaction rate when both par- 
ticles are at the same position, j — i, is then 



\<t>i 



x\nnvu\ A \B r 



h d 



(9) 



This corresponds to the choice of bimolecular reaction 
rate k = A \ B Tb | in the standard RDME. With this choice, 
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when Tb/h —> the RDME may be interpreted as an 
asymptotic approximation of the CRDME. This approx- 
imation is illustrated in Fig. 2. 

In three-dimensions, d — 3, the macroscopic bimolec- 
ular reaction rate k — X\B rb \ also arises as the lead- 
ing order asymptotic expansion as — > 0, A — > 0, or 
D = D A + D B — >• oo of the diffusion limited bimolecu- 
lar reaction rate for the Doi model ([7]), /cooi- In [H] the 
latter was found to be 



■-Doi 



AnDrh ( 1 -^Y? tanh rb V^ 



(Note, the more well-known Smoluchowski diffusion lim- 
ited reaction rate [23 EZ], &Smoi = 4-7rZ3rb, is recov- 
ered in the limit A 4 oo.) As r^^JX/D — > 0, fcDoi ~ 



A(47rrg/3) 



X \B r 



We thus have that the CRDME 



recovers this well- mixed reaction rate as r^/h — > 0. The 
smaller r^^JX/ D, the better the RDME should approxi- 
mate the CRDME for fixed r h /h. 

When modeling three-dimensional biological systems, 
if ?*b\/A/-D is sufficiently small h may simply be chosen 
to accurately model molecular diffusion by a continuous- 
time random walk. In this case, if r^/h <C 1 we may 
approximate the CRDME by the standard RDME. When 
these assumptions break down we need to decrease h and 
incorporate reactions between molecules in neighboring 
voxels with reactive transitions rates X<j>ij. The reaction 
and diffusion processes could potentially be decoupled by 
choosing separate meshes for each (as in [H]). That said, 
this process should be done so as to provide a convergent 
approximation of Q. 
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FIG. 4: Comparison of area fraction when both 
molecules are in the same square, <pu, with area fraction 
when the two molecules are in neighboring squares, 
We see that the area fraction decreases faster as 
r^/h — > when the two molecules are in different 
squares than when they are in the same square. 
Moreover, as r^/h — ¥ the area fraction 4>a approaches 
nrl/h 2 as derived in 



VI. CONCLUSION 

By discretizing the stochastic reaction-diffusion model 
of Doi [6l [7] we have derived a new convergent reaction- 
diffusion master equation. We illustrated our discretiza- 
tion procedure, and the convergence of the survival time 
distribution and mean reaction time in the CRDME, 
for two molecules that undergo the annihilation reaction 
A + B — > 0. While this special case is simplified com- 
pared to realistic biological networks, it should be noted 
that the same reaction rates, X<pij, are derived when 
applying our discretization procedure to the more gen- 
eral multi-particle Doi model 0. Moreover, the same 
type of finite- volume approach we have taken to derive a 
CRDME approximating the Doi model could be applied 
to derive rates for approximating Smoluchowski models. 

We do not show here that the CRDME provides a con- 
vergent discretization in (i-dimensions. We hope to report 
on the convergence in three-dimensions for more realistic 
biological systems in the future. 

Finally, it should be noted that while we have provided 
a method for deriving a convergent reaction-diffusion 
master equation as h — > 0, non-convergent reactive tran- 
sition rates as derived in [T21[]31[5T] may provide more ac- 
curacy in resolving specific statistics of underlying contin- 
uum stochastic reaction-diffusion models for coarse mesh 
widths {h > rb). In particular, several of these works de- 
rive non-convergent reactive transition rates for use in 
the RDME that approximate statistics of Smoluchowski 
diffusion limited reaction models (as opposed to the Doi 
model we have studied). 
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SUPPLEMENTARY MATERIAL 

VII. CALCULATION OF REACTION 
TRANSITION RATES 

Denote by Vi the d-dimensional coordinate axis aligned 
hypercube with sides of length h centered at ih. With 
this definition we may then write Vij = Vi x Vj. We use 
Vi to denote this hypercube in the special case that h = 1. 
Finally, let B rh (x) be the d-dimensional hypersphere of 
radius about x. A convenient representation for 
we subsequently use is 



l — f f l R {\x-y\)dydx (10) 
Tl I ^(aOnVil dx 



Vi 



B rb /h(x)nVj-i dx. 



(11) 



(Here denotes the origin voxel.) Hence we may inter- 
pret <pij as the integral over the center of a hypersphere of 
the volume of intersection between the hypersphere and 



a hypercube. The final equation (11) shows that (f>ij de- 



pends on only two quantities; the separation vector j — i 
and rh/h. Also note that (f>ij will be zero once the sepa- 
ration between all points in voxels i and j is more than 
rt,. As such, in practice it is only necessary to calculate 
4>oj for a small number of voxels about the origin. 

It is desirable to calculate tfiij to near machine precision 
to avoid the introduction of error from the use of incorrect 
reactive transition rates. While this may seem an easy 
task, simply calculating the hypervolume of intersection 
of 7Z and Vij, it should be noted that these are four- 
dimensional (six-dimensional) sets when the molecules 
are in two-dimensions (three-dimensions). Evaluating 



4>ij by directly applying quadrature to ( 10 1 is compli- 
cated by the discontinuous integrand. We have found 
that several standard cubature HI] and Monte Carlo 
methods [19] have difficultly evaluating such integrals in 
reasonable amounts of computing time to high numerical 
precision (absolute errors of near 10~ 12 ). Since the inte- 
gral ( 11 ) has a continuous integrand, which only requires 



the intersection of two-dimensional (three-dimensional) 
sets when the particles are each in two-dimensions (three- 
dimensions), we focus on evaluating <fiij through this rep- 
resentation. 



To evaluate (111 both efficiently and accurately it is 



necessary to calculate the hypervolume of intersection 



given by the integrand, Vj(x) = B rh / h (x) n Vj 



Our 



approach is based on writing this hypervolume as an inte- 
gral and then converting to a boundary integral through 
the use of the divergence theorem. That is, 



Vj(x) 



V ydy, 



y -v(y) dS(y), 



9S rb/h (x) 
1 



+ - / (yr,(y))l B (x) {y)dS{y). (12) 

« JdVj 

Here dM is used to denote the boundary of a manifold M, 
T](y) the outward normal to the boundary hypersurface 
at y, and dS(y) the hypersurface measure at y. 

For simplicity, in the remainder we assume d = 2. In 
this case we have developed a fast method, requiring only 
a few minutes on a modern laptop, that is able to evalu- 



ate ( 11 ) to near machine precision. Vj(x) is evaluated by 



calculating the intersection points of the circle dB rh {x) 
with the square Vj numerically. Once these points are 



known the line integrals in ( 12 ) can be reduced to sums 



of integrals over sub-arcs where the indicator function 
is identically one or zero. These integrals can be evalu- 
ated analytically. Standard adaptive numerical quadra- 
ture methods, such as the dblquad routine in MATLAB, 
are then able to effectively integrate the area of intersec- 
tion function Vj(x). This method was used to generate 
the area fractions in Fig. 4 and all SSA simulations. 
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